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Abstract — The  bidomain  equations  represent  the  most  com¬ 
plete  description  of  cardiac  electrical  activity.  However,  the 
equations  prove  computationally  burdensome  as  the  resulting 
system  of  equations  has  two  entries  per  spatial  node.  This 
paper  examines  the  computational  performance  obtained  by 
decoupling  the  bidomain  equations  into  two  separate  systems 
of  equations,  an  elliptic  equation  for  the  extracellular  poten¬ 
tial,  and  a  parabolic  equation  for  the  transmembraiiie  voltage. 
Each  set  of  equations  was  solved  on  different  grids  with  differ¬ 
ent  time  steps.  For  the  elliptic  problem,  the  performances  of 
direct  and  iterative  solvers  were  bomphred.  For  the  parabolic 
equation,  the  interconnected  cable  method  (ICCM)  was  com¬ 
pared  to  the  finite  element  method  (FEM).  Results  were 
obtained  by  simulating  activity  in  a  3D  slab  of  cardiac  tis¬ 
sue  whose  ionic  currents  were  described  by  modified  Beeler- 
Reuter  equations.  It  Was  found  that  'the  elliptic  equation 
solution  dominated  the  calculation.  Reducing  the  frequency 
of  solution  and/or  halving  the  spatial  resolution  resulted  in 
considerable  speed  tip  while  maintaining  a  reasonable  error. 
Direct  solvers  were  faster  by  a  factor  of  2— 3  and  the  ICCM 
was  about  twice  as  fast  in  solving  the  parabolic  equation  as 
compared  to  the  FEM.  Both  the  elliptic  and  parabolic  equa¬ 
tions  scaled  linearly  with  the  number  of  nodes. 
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a  time  step  on  the  order  of  microseconds  while  the  window 
of  observation  to  determine  the  outcome  of  a  shock  is  on 
the  order  of  hundreds  of  milliseconds.  Thus,  both  spatial 
and  temporal  considerations  contribute  to  the  size  of  the 
problem. 

Many  techniques  are  available  to  solve  the  reaction- 
diffusion  equations  describing  cardiac  electrical  activity. 
The  Interconnected  Cable  Method  (ICCM)  [4],  [5]  is  a  com¬ 
putationally  efficient  method  that  has  been  used  in  mon¬ 
odomain  simulations  of  three  dimensional  cardiac  tissue 
with  fiber  rotation[6].  It  is  based  on  decomposing  the  tissue 
into  a  set  of  cables  which  may  follow  arbitrary  trajectories, 
but,by  itself,  is  not  suitable  for  solving  the  bidomain  equa¬ 
tions.  The  finite  element  method  (FEM)  allows  modeling 
of  complex  geometry  and  has  been  used  to  solve  bidomain 
problems  on  a  whole  rabbit  heart  [7],  but  is  more  computa¬ 
tionally  demanding  than  the  ICCM.  * 

This  study  examines  several  techniques  to  increase  the 
computational  efficiency  of  solving  the  bidomain  equations. 
Benefits  to  be  gained  from  recasting  the  bidomain  equations 


I.  Introduction 

LECTRICAL  shocks  are  the  only  known  therapy  for 
hearts  in  fibrillation.  Without  intervention,  death  will 
quickly  result.  The  mechanisms  underlying  electrical  defib¬ 
rillation,  however,  still  remain  elusive.  Electrical  mapping 
of  the  cardiac  electrical  activity  during  the  delivery  of  de- 
fibrillatory  shocks  is  hindered  by  the  high  shock  strength, 


into  decoupled  elliptic  and  hyperbolic  problems  are  exam¬ 
ined.  With  the  problems  isolated,  each  is  solved  with  dif¬ 
ferent  time  steps  and  on  different  meshes  in  order  to  reduce 
computational  demand.  By  comparing  results  with  the  fully 
coupled  bidomain  solution,  the  simulation  parameters  which 
maximize  computational  speed  while  maintaining  sufficient 
accuracy  are  determined. 

II.  Methods 


while  optical  mapping  techniques  are  limited  to  surface  mea¬ 
surements.  Modeling  is  therefore  needed  to  help  resolve  the 
events  during  defibrillation  in  the  3D  volume  of  the  heart. 

The  bidomain  representation  of  cardiac  tissue  is  the  most 
complete  description  of  cardiac  electrical  activity[l].  It  de¬ 
scribes  both  the  intracellular  and  extracellular  potential 
fields,  linking  them  through  membrane  behavior.  It  has  pre¬ 
dicted  the  appearance  of  shock-induced  virtual  elec.trodes[2] 
which  were  later  confirmed  experimentally  [3].  The  bidomain 
equations  are  computationally  expensive  as  they  require  two 
unknowns  for  each  spatial  node,  resulting  in  large  matri¬ 
ces  which  consume  much  memory  and  require  long  solution 
times. 

Besides  this  inherent  computational  burden  in  solving  the 
bidomain  equations,  simulating  defibrillation  in  the  heart  is 
also  an  intrinsically  large  problem.  A  piece  of  tissue  large 
enough  to  support  fibrillation  must  be  modeled,  on  the  order 
of  centimeters,  and  because  of  the  length  constants  involved, 
it  must  be  discretized  on  the  order  of  hundreds  of  microm¬ 
eters.  Furthermore,  the  kinetics  of  the  sodium  gate  impose 


A.  Governing  Equations 

The  basic  bidomain  equations[l]  relate  the  intracellular 
potential,  fa,  to  the  extracellular  potential,  fa,  through  the 
transmembrane  current  density.  Im : 


fa  —  pljn  Itrans 
V  *  deVrig*  =  film  f*  Arons  A 

J  -C  9Vm  I  / 

±m  '~ym  ^  i  1xon 


(1) 

(2) 

(3) 


where  fa  and  fa  are  respectively  the  intracellular  and  extra¬ 
cellular  conductivity  tensors,  0  is  the  surface  to  volume  ra¬ 
tio  of  the  cardiac  cells,  I(ron5  is  the  transmembrane  current 
stimulus,  Ie  is  an  extracellular  current  stimulus,  Cm  is  the 
capacitance  per  unit,  area,  Vm  is  the  transmembrane  volt¬ 
age  which  is  defined  as  fa -fa,  and  Aon  is  the  current  den¬ 
sity  flowing  through  the  ionic  channels.  The  Beeler-Reuter 
Drouhard-Roberge  model  modified  to  handle  large  voltages 
was  used  as  the  ionic  model[8]  in  this  study.  This  formula¬ 
tion  will  be  referred  to  as  the  coupled  set  of  equations  since 
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the  intracellular  and  extracellular  potentials  are  solved  si¬ 
multaneously. 

By  adding  Eqn.  1  and  Eqn.  2  and  using  the  definition  of 
Vm,  the  equations  can  be  cast  in  a  slightly  different  form 
with  Vm  and  <pe  as  the  independent  variables. 

V  •  ((fi  +  de)V<f>e  =  -  V  -  <nWm  -  Ie  (4) 

V  •  VVjvj  =  \J  .  tTj  V (pc  -f-  f3Im  —  Itrans  (5) 

Equations  4  and  5  are  decoupled  and  solved  sequentially 
as  an  elliptic  problem  (Eqn.  4)  and  a  parabolic  problem 
(Eqn.  5).  Since  the  two  systems  are  now  independent  of 
each  other,  they  can  be  solved  with  different  time  steps  and 
at  different  spatial  resolutions.  This  formulation  will  be 
referred  to  as  the  decoupled  system. 

If  either  the  extracellular  electric  field  can  be  ignored,  or, 
the  ratios  of  the  longitudinal  to  transverse  conductivities  in 
the  intracellular  and  extracellular  domains  are  equal,  the 
bidomain  equations  can  be  replaced  with  the  monodomain 
equation.  This  is  Eqn.  5  with  the  intracellular  conductivity 
tensor  replaced  by  the  the  monodomain  conductivity  ten¬ 
sor,  am-.  which  is  a  function  of  the  bidomain  conductivity 
tensors[9],  pm  =<fi(d:i.  +  £fe)~1<fe. 

B.  Solution  Methods 

To  solve  the  fully  coupled  system  (Eqn.’s  1  and  2)  the 
FEM  was  used  based  on  a  Galerkin  formulation.  In  matrix 
notation  with  a  time  step  of  At,  the  resultant  discretized 
system  is  given  by 

Ke  —  «M  kM  1  \<t>\+l]  _ 
kM  K,  -  kMJ  [^‘+iJ  “ 


where  k  =  j3Cm/At,  M  is  the  FEM  lumped  mass  matrix 
and  K  is  the  FEM  stiffness  matrix.  Both  matrices  were 
computed  using  linear  tetrahedral  elements  with  the  sub¬ 
script  on  K  denoting  whether  the  matrix  was  created using 
Oi (i)  or  <fe(e).  Superscripts  refer  to  the  time  step. 

To  solve  the  elliptic  equation  of  the.  decoupled  system 
(Eqn.  4),  an  FEM  approach  was  also  used: 

K,+e^e  =  Kjv^,  -  Mi*  (7) 

where  the  subscript  of  the  stiffness  matrix  denotes  that  the 
sum  of  the  conductivity  tensors  was  used. 

To  solve  the  parabolic  problem  of  the  decoupled  system 
(Eqn.  5),  two  different  schemes  were  used,  ICCM  and  FEM. 
The  ICCM  solution  utilized  a  semi-implicit  time  integration 
on  a  one-dimensional  linear  grid  to  solve  the  particular  so¬ 
lution^].  The  FEM  solution  utilized  a  forward  Euler  time 
integration: 

vm 1  ~  vm  +  ~(m~1Ki(v^1  —  <f>*)  -Mon  +  i lJ)  (8) 


ICCM  Fine  FEM  Coarse  FEM 


Fig.  1.  .Grids  used  for  different  solution  methods.  AH  grids  represent 
the  same  space.  They  are  shown  in  2D  for  clarity. 


C.  Grid  Generation 

Grid  generation  began  by  constructing  an  ICCM  grid  as 
has  been  described  previously [5]  (see  Figure  1).  Cables  were 
laid  in  sheets  spaced  100  pm  apart.  Within  each  sheet,  ca¬ 
bles  were  parallel  and  all  ran  in  the  y  direction.  The  domain 
spanned  by  the  cables  was  a  right-angled  hexahedron.  Ca¬ 
bles  were  then  discretized  into  100  yzm  long  segments  with 
each  segment  connected  to  a  segment  in  a  neighboring  ca¬ 
ble  through  a  gap  junction  that  was  represented  by  a  fixed 
resistance.  '. 

Two  different  FEM  meshes  were  constructed  from  the  ca¬ 
bles,  a  fine  mesh  and  a  coarse  mesh.  To  construct  the  fine 
three-dimensional  FEM  mesh,  the  centers  of  the  cable  seg¬ 
ments  were  used  as  nodes  from  which  to  construct  first-order 
tetrahedrons.  To  construct  the  coarse  FEM  mesh,  every 
other  point  in  every  other  cable  of  every  other  layer  was 
used  resulting  in  a  mesh  that  was  approximately  one-eighth 
the  size  of  the  finer  mesh.  The  conductivity  tensors  in  the 
tetrahedral  elements  were  defined  by  the  cable  directions  in 
the  ICCM  model.  .  , 

The  effect  of  using  a  coarse  mesh  for  the  elliptic  equa¬ 
tion  and  a  fine  mesh  for  the  parabolic  equation  was  tested. 
Once  a  solution  was  obtained  for  4>e  on  the  coarse  mesh, 
interpolation  was  used  to  assign,  potential  values  at  the  fine 
mesh  nodes  which  were  used  in  the  solution  of  the  parabolic 
equation.  . .  r  ' 

D.  Matrix  Solvers  ■ 

To  solve  the  coupled  set  of  equations  and  the  elliptic  equa¬ 
tion,  both  iterative  and  direct  methods  were  used.  While 
direct  methods  are  generally  faster  when  repeatedly  solv¬ 
ing  the  same  system  of  equations,  they  require  much  more 
memory  since  performing  a  decomposition  on  a  sparse  ma¬ 
trix  preserves  matrix  bandwidth  but  fills  in  the  zero  entries 
between  bands[10].  One  must  therefore  use  iterative  meth- 
ods  on  large  problems  where  performing  a  matrix  decompo¬ 
sition  would  exceed  computer  memory.  The  direct  method 
used  here  was  an  SGI  (Mountain  View,  CA)  supplied  LDLT 
decomposition,  where  L  is  a  lower  triangular  matrix  and  D 
is  a  diagonal  matrix.  After  the  decomposition,  the  system 
was  solved  by  forward  and  backward  substitutions.  A  cus- 


tom  coded  conjugate  gradient  method  with  an  Incomplete 
Cholesky  decomposition  preconditioner  was  used  as- the  it¬ 
erative  solverfll]. 

All  simulations  were  performed  on  an  SGI  Origin  2100 
computer  which  had  350  MHz  MIPS  R12000  processors  and 
4  gigabytes  of  memory.  Times  given  for  simulations  are  CPU 
times  for  a  single  processor. 

III.  Results 

A.  Solution  Methods 

The  CPU  time  taken  to  simulate  25  ms  of  activity  in 
a  rectangular  3D  block  of  cardiac  tissue  with  zero  flux 
boundary  conditions  w'as  measured.  The  block  measured 
1.6  x  0.6  x  0.11  (x  x  y  x  z)  cm  and  was  composed  of  108,031 
intracellular  nodes.  Activity  was  initiated  in  one  corner  of 
the  block  by  applying  a  2  ms  suprathreshold  transmembrane 
stimulus.  The  activity  propagated  out  with  an  ellipsoidal 
wavefront  which  reached  the  x  and  y  edges  at  approximately 
25  ms.  The  CPU  time  to  perform  the  simulation  is  given  in 
Table  I  for  the  various  methods.  AtE/Atp  refers  to  the  ratio 
of  the  time  step  used  for  the  elliptic  equation  solve  over  the 
time  step  used  for  the  parabolic  equation  solve.  Parabolic 
solves  were  always  performed  with  a  time  step  of  10  fis. 
Thus,  0e  was  not  necessarily  updated  as  frequently  as  vra. 
The  value  of  oo  corresponds  to  solving  the  monodomain 
equation.  Fine  and  Coarse  refer  to  the  discretization  of  the 
FEM  grid  on  which  the  elliptic  equation  was  solved.  Itera¬ 
tive  refers  to  solving  the  elliptic  equation  on  a  fine  grid  using 
the  conjugate  gradient  method. 

TABLE  I 

CPU  TIME  TO  SIMULATE  25  MS  FOR  VARIOUS  METHODS. 


Method 

AtE/Atp 

Fine 

Coarse 

Iterative 

ICCM 

1 

3989.67 

941.20 

8517.56 

2 

2270.54 

796.63 

5909.88 

4 

1396.46 

629.02 

4115.56 

10 

878.25 

602,41 

2761.32 

00 

398.42 

- 

'  - 

FEM 

1 

4868.44 

1171.68 

10071.34 

2 

2515.26 

979.92 

7111.40 

4 

1643.83 

884.71 

5018.22 

10 

1123.82 

819.91 

3404.27 

00 

757.88 

-  - 

- 

Coupled 

13489.57 

— 

65867.97 

Solving  the  coupled  equation  was  much  slower  than  solv¬ 
ing  the  decoupled  system.  The  ICCM  method  was  ap¬ 
proximately  twice  as  fast,  as  the  FEM  method  in  solving 
the  parabolic  problem.  An  iterative  solver  was  slightly 
more  than  twice  as  slow  to  solve  the  elliptic  problem  at 
a  AtE/Atp  of  one  and  became  three  times  as  slow  when 
AtE /Atp  was  ten.  This  was  due  to  the  larger  change  in 
4>e  when  computed  less  often,  thereby  requiring  more  itera¬ 
tions  for  convergence.  Increasing  AtE/Atp  from  one  to  ten 
increased  the  average  number  of  iterations  per  elliptic  solve 
from  52  to  94. 


Fig.  2.  Time  spent  in  the  parabolic  solution  for  a  25  ins  long  sim¬ 
ulation  as  a  function  of  problem  si2e.  The  size  of  the  problem 
was  increased  by  increasing  the  number  of  layers  (#layers)  or  by 
increasing  the  size  of  each  layer  (layer  size). 

The  elliptic  solve  was  the  most  costly  part  of  the  solution. 
Reducing  the  number  of  times  that  the  elliptic  equation  is 
solved  (increasing  AtE / Atp)  or  reducing  the  size  of  the  el¬ 
liptic  problem  (using  a  coarse  grid)  greatly  decreased  the 
simulation  time.  Monodomain  solutions  were  obviously  the 
quickest. 

B.  Problem  Size 

The  dependence  of  simulation  time  on  problem  size  was 
next  ascertained.  The  problem  size  was  varied  in  two  dif¬ 
ferent  manners.  The  first  manner  was  to  increase  the  z  di¬ 
mension  which  simply  increased  the  number  of  layers  in  the 
block.  This  resulted  in  an  increase  in  the  number  of  nodes 
while  preserving  the  bandwidth  of  any  matrices.  The  second 
method  involved  scaling  the  size  of  each  layer  in  the  x  and 
y  directions  by  the  same  factor  while  keeping  the  number 
of  layers  constant.  This  method  resulted  in  matrices  whose 
bandwidth  increased  with  the  scaling  factor.  The  effect  of 
problem  size  was  broken  down  into  the  effect  on  the  time  of 
the  parabolic  solution  (Fig.  2)  and  the  effect  on  the  time  of 
the  elliptical  solution  for  a  coarse  gird  (Fig.  3). 

The  ICCM  method  was  about  twice  as  fast  as  the  FEM 
method  in  solving  the  parabolic  equation.  The  ICCM 
method  scaled  linearly  with  problem  size,  regardless  of  how 
the  problem  size  was  increased.  The  curves  for  the  ICCM 
method  were  on  top  of  each  other  since  bandwidth  is  not  a 
consideration  for  this  method. 

The  FEM  showed  a  linear  increase  in  CPU  time  with 
problem  size.  At  the  smaller  sizes,  the  manner  in  which  the 
problem  wras  increased  in  size  did  not  affect  the  computation 
time.  However,  for  very  large  problems,  more  than  200,000 
nodes,  increasing  the  layer  size  caused  a  larger  increase  in 
computation  time  compared  to  an  increase  in  the  number 
of  layers. 

The  effect  of  increasing  the  problem  size  for  the  ellipti¬ 
cal  problem  was  similar  to  the  one  obtained  for  the  FEM 
parabolic  solve.  The  manner  in  which  the  size  was  increased 


Fig.  3.  Time  spent  solving  the  elliptic  equation  in  a  simulation  of  25  as 
a  function  of  problem  size.  The  number  of  nodes  in  the  model  was 
increased  by  either  increasing  the  number  of  layers  in  the  model 
or  by  increasing  the  size  of  each  layer. 

was  only  significant,  for  very'  large  problems. 

IV.  Discussion 

This  paper  analyzed  the  time  to  compute  a  simulation  us¬ 
ing  the  bidomain  model  for  various  methods  and  techniques. 
Solving  the  decoupled  System  of  the  bidomain  equations  is 
computationally  advantageous.  By  splitting  the  problem 
into  two,  each  part  can  be  solved  independently.  This  split¬ 
ting  leads  to  two  smaller  problems  whose  total  work  is  less 
than  that  required  to  solve  the  coupled  system.  Further¬ 
more..  there  is  a  large  savings  in  memory'  since  the  matrices 
to  be  constructed  are  much  smaller. 

In  splitting  the  bidomain  equations  into  two  parts,  it  was 
found  that  the  elliptic  equation  was  much  more  costly  to 
solve  than  the  parabolic,  equation.  Thus,  computational 
speed  was  increased  by  a  combination  of  solving  the  elliptic 
problem  on  a  grid  with  coarser  discretization  or  by  solving 
this  equation  less  often  than  the  parabolic  equation.  Us¬ 
ing  a  coarse  grid  had  a  similar  speedup  to  decreasing  the 
frequency  at  which  the  elliptic  equation  was  solved  to  one 
tenth.  Performing  the  elliptic  solve  at  less  than  this  fre¬ 
quency  lead  to  large  errors  in  computation.  Combining  a 
coarse  grid  with  periodic  elliptical  solves  only'  increased  the 
performance  marginally  beyond  a  AtE/Atp  value  of  4. 

The  ICCM  was  found  to  be  about  twice  as  fast  as  the 
FEM  for  solving  the  parabolic  portion  of  the  bidomain  prob¬ 
lem.  However,  when  the  elliptic  portion  and  ionic,  compu¬ 
tation  are  factored  into  the  problem,  the  savings  in  compu¬ 
tation  time  between  using  the  FEM  and  ICCM  reduces  to 
about  30%.  Also,  this  savings  of  30%  will  diminish  if  a  more 
detailed  ionic  model  is  used.  The  ionic  model  used  here  was 
very  simple,  only  having  seven  state  variables,  while  recent 
models  have  over  20[12],  The  added  complexity  in  setting 
up  the  problem  to  utilize  both  an  ICCM  grid  and  FEM  grid 
may  not  be  justified  for  an  increase  in  performance  of  only 
10%. 


The  direct  solver  for  the  elliptic  elliptic  was  2-3  times 
faster  than  the  iterative  solver.  It  is  therefore  only  reason¬ 
able  to  use  the  iterative  solver  if  memory  is  an  issue  and 
the  decomposed  matrix  does  not  fit  into  memory'.  A  decom¬ 
posed  sparse  matrix  can  have  10-20  times  as  many  nonzero 
entries  as  the  original  matrix  due  to  fill  in.  This  an  addi¬ 
tional  reason  to  solve  the  elliptic  equation  on  a  coarse  FEM 
grid  since  memory  for  an  LDl?  decomposition  will  be  re¬ 
duced  by  a  factor  of  eight. 

V.  Conclusions 

Decoupling  the  bidomain  equations  into  an  elliptic  and 
parabolic  equation  offer  computational  advantages  which 
can  be  exploited  to  solve  much  larger  problems  in  much  less 
time.  The  elliptic,  problem  was  computationally  more  ex¬ 
pensive  as  well  as  requiring  more  memory.  CPU  times  were 
greatly  reduced  by  solving  it  on  a  coarser  spatial  grid  and  at 
fewer  instances  in  time  while  keeping  errors  within  reason¬ 
able  bounds.  For  solving  the  parabolic  problem,  the  ICCM 
was  approximately  twice  as  fast  as  the  FEM.  Due  to  the 
sparsity'  of  the  problem,  both  the  elliptic,  and  parabolic  prob¬ 
lems  scaled  linearly  with  problem  size.  Finally,  for  problems 
which  fit  into  memory,  direct  methods  are  two  to  three  times 
faster  than  iterative  methods. 
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